import numpy as np
from scipy.optimize import minimize
import pandas as pd
from IPython.display import clear_output
from scipy.optimize import fsolve
import time
import matplotlib.pyplot as plt
import networkx as nx
import copy

from InitProblem_library import RKM_Dyn_System



def Delta_Kron(i, j):
    
    """
    The Kronecker delta function
    """
    
    if i == j:
        return 1
    else:
        return 0

def F_y_i_j(a, f, i, j, y, u, lambda_, Arguments):
    
    """
    This function calculate the partial derivative of the funtion F_{a, f} with repsect to the variable y_{i, j} 
    at the fixed time moment (see formula ... in the main manuscript)
    
    Input:
    a         - the first index of F_{a, f} (opinion dimension)
    f         - the second index of F_{a, f} (type dimension)
    i         - the first index of y_{i, j} (opinion dimension)
    j         - the second index of y_{i, j} (type dimension)
    y         - the matrix of state variables
    u         - the matrix of control variables
    lambda_   - the matrix of adjoint variables
    Arguments - hyperparameters
    """    
    
    m = Arguments['m']
    M = Arguments['M']
    
    n = Arguments['n']  # the fraction of ordinary agents
    
    Transition_Matrices = Arguments['Transition_Matrices']
    
    w = Arguments['w']  # the opinion weights (integral)
    
    #--------------------------------------------------------
    
    Derivative = 0  # the output variable
    
    n_f = u[:, f].sum()  # the fraction of stubborn agents that communicate with type f agents at the moment
            
    for l in range(m):
                
        for r in range(M):
                    
            Derivative = Derivative + Delta_Kron(j, f) * y[l, r] * Transition_Matrices[f'{j}-{r}'][i][l, a]
                    
        Derivative = Derivative + Delta_Kron(j, f) * u[l, j] * Transition_Matrices[f'{j}-{M}'][i][l, a]
                    
    for s in range(m):
                    
        Derivative = Derivative + y[s, f] * Transition_Matrices[f'{f}-{j}'][s][i, a]
                    
    Derivative = Derivative - Delta_Kron(i, a) * Delta_Kron(j, f) * (n + n_f)
    
    return Derivative    
    

def G_i_j(i, j, y, u, lambda_, Arguments):
    
    """
    This function calculate the component G_{i, j} of the velocity matrix G (see formula (...) in the 
    main manuscript) multiplied by (-1).
    
    Input:
    i         - the first index of F_{a, f} (opinion dimension)
    j         - the second index of F_{a, f} (type dimension)
    y         - the matrix of state variables
    u         - the matrix of control variables
    lambda_   - the matrix of adjoint variables
    Arguments - hyperparameters
    """
    
    m = Arguments['m']
    M = Arguments['M']
    
    n = Arguments['n']  # the fraction of ordinary agents
    
    Transition_Matrices = Arguments['Transition_Matrices']
    
    #--------------------------------------------------------
    
    Velocity = 0  # the output variable
    
    w = Arguments['w']  # the opinion weights (integral)
    
    for a in range(m):
        for f in range(M):
            
            Velocity = Velocity - lambda_[a, f] * F_y_i_j(a, f, i, j, y, u, lambda_, Arguments)
    
    Velocity = Velocity + w[i]
    
    return (-1) * Velocity



def G(y, u, lambda_, Arguments):
    
    """
    This function calculate the velocity matrix G (see formula (...) in the 
    main manuscript) multiplied by (-1).
    
    Input:
    y         - the matrix of state variables
    u         - the matrix of control variables
    Arguments - hyperparameters
    """
    
    m = Arguments['m']
    M = Arguments['M']
    
    #--------------------------------------------------------
    
    Velocities = np.zeros((m, M))  # the output variable
    
    for i in range(m):
        for j in range(M):
            Velocities[i, j] = G_i_j(i, j, y, u, lambda_, Arguments)
            
    return Velocities



def RKM_Iteration_Euler_Lagrange_System(y_previous, u_previous, Lambda_previous, Arguments):
    
    """
    This function performs one Runge-Kutte iteration for system (...)
    
    Input:
    y_previous - the values of the state matrix at the previous time moment
    u_previous - the values of the control matrix at the previous time moment
    Arguments  - hyperparameters
    
    Output:    - the values of the state matrix at the next time moment
    """
    
    m = Arguments['m']
    M = Arguments['M']
    
    # The elements of the Butcher tableau
    
    #b = Arguments['b']
    #s = Arguments['s']
    
    # The grid
    
    Step = Arguments['Step']
    
    #--------------------------------------------------------
        
    K1 = Step * G(y_previous, u_previous, Lambda_previous, Arguments)
    K2 = Step * G(y_previous, u_previous, Lambda_previous+K1/2, Arguments)
    K3 = Step * G(y_previous, u_previous, Lambda_previous+K2/2, Arguments)
    K4 = Step * G(y_previous, u_previous, Lambda_previous+K3, Arguments)
        
    Lambda_new = Lambda_previous + K1/6 + K2/3 + K3/3 + K4/6
 
    return Lambda_new

    

def RKM_Euler_Lagrange_System(y_t, u_t, Arguments, iteration_counter):
    
    """
    This function realizes the Runge-Kutta procedure for the Cauchy problem (...). 
    
    Note: instead of solving problem (...), this function concetrates on (...)
    
    Input:
    y_t               - the array of state matrices (y (\tau) in the main manuscript)
    u_t               - the array of control matrices (u (\tau) from the main manuscript) 
    Arguments         - hyperparameters
    iteration_counter - the number of the current FBS method iteration
    
    Output:
    lambda_t          - the array of adjoint matrices (\lambda (\tau) in the main manuscript)
    """
    
    m = Arguments['m']
    M = Arguments['M']
    
    # Initial point (defined at the right edge of the interval but ...)
    
    v = Arguments['v']  # the opinion weights (terminal)
        
    # The grid
    
    Step = Arguments['Step']
    Grid = Arguments['Grid']
        
    #------------------------------------------------------------------
    
    # Initialize the value of the adjoint function on the right end of the interval [ \tau_{0}, \tau_{1} ]
    
    lambda_0 = np.ones((m, M))
    
    for i in range(m):
        
        lambda_0[i, :] = (-1) * v[i]
        
    # Initialize the output variable
    
    lambda_t = []
    lambda_t.append(lambda_0)    
    
    # The number of Runge-Kutta iterations
    
    N_Iter = len(Grid)-1
    
    # Start iterations
    
    for i in range(N_Iter):
        
        # Display the current iteration
        
        #print(f'Iteration {iteration_counter} ... Solving Euler-Lagrange system: {i}/{N_Iter}')
        #clear_output(wait=True)
        
        # We will solve a reversed system (in other words, we proceed the system from thr right end point to 
        # the left one)
        
        y_previous = y_t[-i]  #!!
        u_previous = u_t[-i]  #!!!!!!!
        
        lambda_previous = lambda_t[i]
        
        lambda_new = RKM_Iteration_Euler_Lagrange_System(y_previous, u_previous, lambda_previous, Arguments)
        
        
            
        # Check of the ajoint function is nonpositive
        
        #if len(np.where(Lambda_new > 0)[0]) > 0:
            
        #    raise ValueError(f'The adjoint function inlcudes positive components!')
            
        #else:
            
        #    pass
        
        lambda_t.append(lambda_new)
    
    # As an output, we give a reversed array (see Appendix B in the main manuscript for details)
        
    lambda_t.reverse()
    
    return lambda_t



#####################################################################################################################
#####################################################################################################################
#####################################################################################################################


def H_u_i_j(i, j, y, lambda_, Arguments):
    
    """
    This function computes the partial derivative of the Hamiltonian function with repsect to the control variable u_{i, j} 
    (see formula (...) in the main manuscript) at the fixed point of the grid
    
    Input:
    i          - the first index of the variable u_{i, j} (opinion dimension)
    j          - the second index of the variable u_{i, j} (type dimension)
    y         - the matrix of state variables
    lambda_   - the matrix of adjoint variables
    Arguments - hyperparameters
    """
    
    m = Arguments['m']
    M = Arguments['M']
    
    Transition_Matrices = Arguments['Transition_Matrices']
    
    #--------------------------------------------------------
    
    Derivative = 0
    
    for a in range(m):
        
        for s in range(m):
            Derivative = Derivative + lambda_[a, j] * y[s, j] * Transition_Matrices[f'{j}-{M}'][s][i, a]
            
        Derivative = Derivative - lambda_[a, j] * y[a, j]
    
    return Derivative



def H_u(y, lambda_, Arguments):
    
    """
    This function computes the partial derivatives of the Hamiltonian with repsect to the control variables u 
    (see formula (...) in the main manuscript) at the fixed point of the grid
    
    Input:
    y         - the matrix of state variables
    lambda_   - the matrix of adjoint variables
    Arguments - hyperparameters
    """
    
    m = Arguments['m']
    M = Arguments['M']
    
    #--------------------------------------------------------
    
    Derivatives = np.zeros((m, M))
    
    for i in range(m):
        for j in range(M):
            Derivatives[i, j] = H_u_i_j(i, j, y, lambda_, Arguments)
    
    return Derivatives


def H_u_i_j_tau(i, j, I_j, y, lambda_, Arguments):
    
    """
    This function computes the coefficients of the formula (...) in the Online Supplementary Materials file
    
    Input:
    i          - the first index of the variable u_{i, j} (opinion dimension)
    j          - the second index of the variable u_{i, j} (type dimension)
    I_j        - the set of indices I_{j} (see ... in the Online Supplementary Materials file)
    y          - the matrix of state variables
    lambda_    - the matrix of adjoint variables
    Arguments  - hyperparameters
    
    Output:
    A_i_j_dict - the dictionary that stores the coefficients
    
    Note: names of coefficients
    
    l_j        - the coefficient A_{l, j}
    free       - the coefficient A_{i, j} (free term)
    """
    
    m = Arguments['m']
    M = Arguments['M']
    
    n = Arguments['n']
    
    w = Arguments['w']
    
    Transition_Matrices = Arguments['Transition_Matrices']
    
    #------------------------------------------------------------------
    
    A_i_j_dict = {}
    
    for l in I_j:
    
        A_l_j = 0

        for a in range(m):
            for z in range(m):

                for k in range(m):

                    A_l_j = A_l_j + (-1)*lambda_[z, j]*Transition_Matrices[f'{j}-{M}'][a][l, z]*y[k, j]*Transition_Matrices[f'{j}-{M}'][k][i, a]

                A_l_j = A_l_j + lambda_[z, j]*Transition_Matrices[f'{j}-{M}'][a][l, z]*y[a, j]

                A_l_j = A_l_j + lambda_[z, j]*Delta_Kron(a, z)
                
                
        for a in range(m):
            for k in range(m):
                for s in range(m):
                        
                    A_l_j = A_l_j + lambda_[a, j]*Transition_Matrices[f'{j}-{M}'][k][i, a]*y[s, j]*Transition_Matrices[f'{j}-{M}'][s][l, k]
                        
            
            for s in range(m):
                    
                A_l_j = A_l_j + (-1)*lambda_[a, j]*y[s, j]*Transition_Matrices[f'{j}-{M}'][s][l, a]
                    
            for a in range(m):
                for k in range(m):
                    
                    A_l_j = A_l_j + (-1)*lambda_[a, j]*y[k, j]*Transition_Matrices[f'{j}-{M}'][k][i, a]
                    
            A_l_j = A_l_j + lambda_[a, j]*y[a, j]
                    
                
        A_i_j_dict[f'{l}-{j}'] = A_l_j
        
    Free_Term = 0
    
    for a in range(m):
        
        factor_1 = 0
        
        factor_1 = factor_1 + w[a]
        
        for z in range(m):
            for f in range(M):
                for l in range(m):
                    for r in range(M):
                        
                        factor_1 = factor_1 + (-1)*lambda_[z, f]*Delta_Kron(j, f)*y[l, r]*Transition_Matrices[f'{f}-{r}'][a][l, z]
                        
                for s in range(m):
                    
                    factor_1 = factor_1 + (-1)*lambda_[z, f]*y[s, f]*Transition_Matrices[f'{f}-{j}'][s][a, z]
                    
                factor_1 = factor_1 + Delta_Kron(a, z)*Delta_Kron(j, f)*n
        
        factor_2 = 0
        
        for k in range(m):
            
            factor_2 = factor_2 + y[k, j]*Transition_Matrices[f'{j}-{M}'][k][i, a]
        
        factor_2 = factor_2 - y[a, j]
        
        Free_Term = Free_Term + factor_1*factor_2
        
    for a in range(m):
        for k in range(m):
            for s in range(m):
                for r in range(M):
                    Free_Term = Free_Term + lambda_[a, j]*Transition_Matrices[f'{j}-{M}'][k][i, a]*y[s, j]*y[l, r]*Transition_Matrices[f'{j}-{r}'][s][l, k]
            
            Free_Term = Free_Term + lambda_[a, j]*Transition_Matrices[f'{j}-{M}'][k][i, a]*(-1)*y[k, j]*n
            
        Free_Term = Free_Term + lambda_[a, j]*y[a, j]*n
            
    A_i_j_dict['free'] = Free_Term     
    
    return A_i_j_dict


def Singular_Control(Max_Derivatives, y, lambda_, Arguments):
    
    """
    This function maximizes the Hamiltonian at the fixed point of the grid in the case of singular control (when there are 
    several partial derivatives \frac { \partial H} { \partial u_{i, j} }) that have the largest value (henceforth - 
    max derivatives)
    
    Input:
    Max_Derivatives - the indices of max derivatives (2D-array, first column - row indices, second - column indices,
                      each row corresponds to one partial derivative)
    y               - the matrix of state variables
    lambda_         - the matrix of adjoint variables
    Arguments       - hyperparameters
    
    Output:
    u               - the obtained estimation of matrix of control variables 
    """
    
    m = Arguments['m']
    M = Arguments['M']
    
    n = Arguments['n']
    
    Max_Derivatives = np.array(Max_Derivatives)
    
    N_Max_Derivatives = Max_Derivatives.shape[0]
    
    #------------------------------------------------------------------
    
    # For each j we compute I_{J}
    
    I_j_dict = {}

    for j in np.unique(Max_Derivatives[:, 1]):

        I_j = np.unique(Max_Derivatives[Max_Derivatives[:, 1] == j][:, 0])

        I_j_dict[f'{j}'] = I_j


    # To obtain a linear system (...) (see main manuscript) we need to go through all of the equations (...)
    # and retrieve the system coefficients (in dictionary form)

    Linear_Sys_Dicts = []

    for k in range(N_Max_Derivatives-1):
        
        # The left part of equation (...)
        
        i_left = Max_Derivatives[k, 0]
        j_left = Max_Derivatives[k, 1]

        A_i_j_dict_left_part = H_u_i_j_tau(i_left, j_left, I_j_dict[f'{j_left}'], y, lambda_, Arguments)
        
        # The right of equation (...)

        i_right = Max_Derivatives[k+1, 0]
        j_right = Max_Derivatives[k+1, 1]

        A_i_j_dict_right_part = H_u_i_j_tau(i_right, j_right, I_j_dict[f'{j_right}'], y, lambda_, Arguments)
        
        # Transfer all coefficients in the left part
        
        for key in A_i_j_dict_right_part.keys():

            if A_i_j_dict_left_part.get(key):

                A_i_j_dict_left_part[key] = A_i_j_dict_left_part[key] - A_i_j_dict_right_part[key]

            else:

                A_i_j_dict_left_part[key] = (-1) * A_i_j_dict_right_part[key]

        Linear_Sys_Dicts.append(A_i_j_dict_left_part)


    # Prepare linear system representation (Ax=b)
    
    A = np.zeros((N_Max_Derivatives, N_Max_Derivatives))

    b = np.zeros(N_Max_Derivatives)

    # Fill A and b
        
    for l in range(len(Linear_Sys_Dicts)):

        Equation_Dict = Linear_Sys_Dicts[l]

        b[l] = Equation_Dict['free']

        for k in range(N_Max_Derivatives):

            i = Max_Derivatives[k, 0]
            j = Max_Derivatives[k, 1]

            if Equation_Dict.get(f'{i}-{j}'):

                A[l, k] = Equation_Dict[f'{i}-{j}']

            else:

                pass
    
    # The last equation in the system reflects the control constraint: 

    A[-1] = 1
    b[-1] = 1-n
    
    # Check if A is singular
    
    if np.linalg.det(A) != 0:
    
        # Solve the system

        solution = np.linalg.solve(A, b) 
        
        # Check the solution
        
        if (solution < 0).sum() > 0:
            
            u = np.zeros((m, M))

            for k in range(N_Max_Derivatives):

                i = Max_Derivatives[k, 0]
                j = Max_Derivatives[k, 1]

                u[i, j] = (1-n)/N_Max_Derivatives               
            
            #raise ValueError(f'Solution of the system for finding control variables includes negative values: {A}, {b}, {solution}!')
            
        else:  # all is OK
            
            u = np.zeros((m, M))

            for k in range(N_Max_Derivatives):

                i = Max_Derivatives[k, 0]
                j = Max_Derivatives[k, 1]

                u[i, j] = solution[k]    
            
    else:
        
        u = np.zeros((m, M))
        
        for k in range(N_Max_Derivatives):
            
            i = Max_Derivatives[k, 0]
            j = Max_Derivatives[k, 1]
            
            u[i, j] = (1 - n) / N_Max_Derivatives
            
        
        #raise ValueError(f'Matrix of the system for finding control variables is singular!')
        
    
    return u



def Maximize_Hamiltonian_Point(y, lambda_, Arguments):
    
    """
    This function maximizes the Hamiltonian at the fixed point of the grid by adjusting the matrix of control variables
    
    Input:
    y         - the matrix of state variables
    lambda_   - the matrix of adjoint variables
    Arguments - hyperparameters
    
    Output:
    u         - the estimated matrix of control variables    
    """
    
    m = Arguments['m']
    M = Arguments['M']
    
    n = Arguments['n']  # the fraction of ordinary agents
    
    #------------------------------------------------------------------
    
    Derivatives = H_u(y, lambda_, Arguments)  # compute all partial derivatives with respect to u
    
    #print(Derivatives)
    
    Max_Value = Derivatives.max()
    
    Max_Derivatives = []
    
    for i in range(m):
        for j in range(M):
            
            if Derivatives[i, j] == Max_Value:
                
                Max_Derivatives.append([i, j])
                
    #Max_Derivatives = np.array(Max_Derivatives)
                
    if len(Max_Derivatives) == 1:  # make use of Proposition ... in the main manuscript
        
        #print(Max_Derivatives[0][0], Max_Derivatives[0][1])
        
        #print('')
        #print('--------------------')
        #print('')
        
        #print(u)
        
        u = np.zeros((m, M))
        
        u[Max_Derivatives[0][0], Max_Derivatives[0][1]] = 1 - n
        
        #print(u)
        
    else:  # singular control - make use of Proposition ... in the main manuscript
        
        u = Singular_Control(Max_Derivatives, y, lambda_, Arguments)
        
        #print(f'Singular control: {Max_Derivatives}!')
        
        #return Max_Derivatives
        
        #raise ValueError(f'Singular control: {Max_Derivatives}!')
    
    #print(u)
    
    #print('')
    
    return u



def Maximize_Hamiltonian(y_t, u_t, lambda_t, Arguments):
    
    """
    This function maximizes the Hamiltonian pointwise by adjusting the control matrix 
    at this point
    
    Input:
    y_t               - the array of state matrices (y (\tau) in the main manuscript)
    u_t               - the array of control matrices (u (\tau) from the main manuscript) 
    lambda_t          - the array of adjoint matrices (\lambda (\tau) in the main manuscript)
    Arguments         - hyperparameters
    
    Output:
    u_t_new           - a new estimation of the control function u (\tau) (an array of control matrices) 
    """
    
    m = Arguments['m']
    M = Arguments['M']
    
    # The grid
    
    Grid = Arguments['Grid']
    
    #------------------------------------------------------------------
    
    # The output variable
    
    u_t_new = []
    
    # Go through the grid points
    
    for i in range(len(Grid)-1):
        
        #print(i)
        
        u_t_new.append(Maximize_Hamiltonian_Point(y_t[i+1], 
                                                  lambda_t[i], 
                                                  Arguments))
    
    """
    # For the last point of the grid, the value of the control function makes no sense
    # we set it to be equal to the previous value
    
    u_t_new.append(u_t_new[-1])
    """
    
    return u_t_new

def Objective_Functional(y_t, Arguments):
    
    """
    For a given array of state matrices (y (\tau) in the main manuscript), 
    this function computes the value of the objective
    functional (see formula (...) in the main manuscript)
    
    Input:
    y_t                  - the array of state matrices (y (\tau) in the main manuscript)
    Arguments            - hyperparameters  
    
    Output:
    
    (Integral, Terminal) - the courtege of the integral and terminal terms of the objective functional
    """
    
    # The grid
    
    Grid = Arguments['Grid']
    
    Step = Arguments['Step']
    
    # Opinion weight vectors
    
    w = Arguments['w']  # integral
    v = Arguments['v']  # terminal
    
    #------------------------------------------------------------------
    
    Integral = 0  
    
    # Compute the integral term
    
    for i in range(len(Grid)-1):
        
        Integral = Integral + ((y_t[i]).sum(axis=1)*w).sum() * Step
        
    # Calculate the terminal term
    
    Terminal = ((y_t[-1]).sum(axis=1)*v).sum()
    
    return (Integral, Terminal)


def FBSM(y_0, u_t_0, Arguments):
    
    """
    This function realizes the FBS method
    
    Input:
    y_0                                               - the starting point of the system (see formula (4) from 
                                                        the main manuscript)
    u_t_0                                             - the initial guess to the array of control matrices (u_{0} (\tau) 
                                                        from the main manuscript) 
    Arguments                                         - hyperparameters
    
    Output:
    (u_t_estimations, Functional_estimations, Status) - the triplet:
    
                                                        (i)   - consecutive estimations of the array of control matrices 
                                                        (ii)  - correpsonding values of the objective functional
                                                        (iii) - convergence status (possible values (see Online 
                                                                                    Supplementary matrials for details):
                                                                                                     (a) - convergence
                                                                                                     (b) - cycle of length 
                                                                                                           k
                                                                                                     (c) - no convergence
                                                                                    )
    """
    
    # Initialize output variables
    
    u_t_estimations = []  # here we will store consecutive algorithm outputs (as lists of 2D arrays)
    u_t_estimations_2D_array = []  # here we will store consecutive algorithm outputs (as 2D arrays)
    Functional_estimations = []  # here we will store the values of the objective functional for each 
                                 # control function estimation
    
    u_t_estimations.append(u_t_0)  # add the initial guess as a first estimation
    u_t_estimations_2D_array.append(np.hstack(u_t_0))
    
    #------------------------------------------------------------------
    
    # start FBSM iterations
    
    iteration_counter = 0
    
    while True:  
        
        iteration_counter = iteration_counter + 1
        
        # Take the previous estimation of the control function
        
        u_t_previous = u_t_estimations[-1]
        
        # Calculate the state function y (\tau)
        
        #y_t = SCARDO_Dyn_System(y_0, u_t_previous, Arguments, iteration_counter)  # solving dynamic system via agent-based approach
        y_t = RKM_Dyn_System(y_0, u_t_previous, Arguments, iteration_counter)
        
        # Calculate the corresponding value of the objective functional
        
        Functional_estimations.append(Objective_Functional(y_t, Arguments))
        
        # Calculate the adjoint function \lambda (\tau)
        
        lambda_t = RKM_Euler_Lagrange_System(y_t, u_t_previous, Arguments, iteration_counter)
    
        # Obtain new estimation to the control function u (\tau)
    
        u_t_new = Maximize_Hamiltonian(y_t, u_t_previous, lambda_t, Arguments)
        
        # Insert new estimation into record arrays
        
        u_t_estimations.append(u_t_new)
        u_t_estimations_2D_array.append(np.hstack(u_t_new))
        
        # Compare new estimation vs previous ones
        
        for i in range(1, len(u_t_estimations)):
            
            """
            print(u_t_estimations_2D_array[-1])
            print('')
            print(u_t_estimations_2D_array[-1-i])
            print('')
            print(u_t_estimations_2D_array[-1] - u_t_estimations_2D_array[-1-i])
            raise ValueError(f'The adjoint function of the output inlcudes positive components!')
            """
            
            Amount_of_Difference = abs(u_t_estimations_2D_array[-1] - u_t_estimations_2D_array[-1-i]).sum()
            
            #N_differences = (u_t_estimations_2D_array[-1] != u_t_estimations_2D_array[-1-i]).sum()
            
            #print(f'Iteration {iteration_counter} ... Check new estimation ... The number of differences with previous estimation: {N_differences}')
            
            #time.sleep(1/4)
            
            if Amount_of_Difference < 0.005:  
                
                # Need to compute the value of the functional for the last estimation
                
                u_t_previous = u_t_estimations[-1]
                
                # Calculate the state function y (\tau)
                
                #y_t = SCARDO_Dyn_System(y_0, u_t_previous, Arguments, iteration_counter)  # solving dynamic system via agent-based approach
                y_t = RKM_Dyn_System(y_0, u_t_previous, Arguments, iteration_counter)
                
                # Calculate the corresponding value of the objective functional
                
                Functional_estimations.append(Objective_Functional(y_t, Arguments))
                
                if i == 1:  # convergence
                    
                    Status = 'Convergence!'
                    
                    clear_output(wait=True)
                    print(f'Complete! Status: {Status}')
                    
                    """
                    for lambda_ in lambda_t:
                    
                        if len(np.where(lambda_ > 0)[0]) > 0:

                            raise ValueError(f'The adjoint function of the output inlcudes positive components!')

                        else:

                            pass
                            
                    """
                    
                    return (u_t_estimations, Functional_estimations, Status)  
                
                else:  # cycle
                    
                    Status = f'Cycle of lenght {i}!'
                    
                    clear_output(wait=True)
                    print(f'Complete! Status: {Status}')
                    
                    return (u_t_estimations, Functional_estimations, Status)  
                    
            else:  # continue the algorithm
                
                if iteration_counter > 300:  # the maximal number of iterations
                    
                    Status = 'No convergence!'
                    
                    clear_output(wait=True)
                    print(f'Complete! Status: {Status}')
                    
                    return (u_t_estimations, Functional_estimations, Status)   
                    
                else:
                    
                    pass

                
                




